Setup

# Remove all from the current workspace
rm(list = ls())

# Set the working directory
setwd('C:\\Users\\srica\\Desktop\\Radio\\projekt\\PDR Tema 9')

# Create a directory named 'plots' if it doesn't already exist
if (!file.exists('plots')) dir.create('plots')

# Create a directory named 'results' if it doesn't already exist
if (!file.exists('results')) dir.create('results')

# Import the 'stringr' library
library(stringr)

Constants

# Set the necessary constants
constants <- data.frame(
  fiPole     = 78.3 * pi / 180,
  lambdaPole = 291 * pi / 180,
  ippHeight   = 350,
  earthRadius = 6378,
  lightSpeed  = 2.99792458e+08
)
print(constants)
##     fiPole lambdaPole ippHeight earthRadius lightSpeed
## 1 1.366593   5.078908       350        6378  299792458

User coordinates

# receiver (user) cooridinates (Pevex Kukuljanovo)
userFiLambda <- data.frame(
  fiUser = 45.331360488178596*(pi/180), # longitude
  lambdaUser = 14.509391791906435*(pi/180) # latitude
)
print(userFiLambda)
##      fiUser lambdaUser
## 1 0.7911815  0.2532367
# real Cartesian receiver (user) coordinates
userCoords <- data.frame(
  X = constants$earthRadius * cos(userFiLambda$lambdaUser) * cos(userFiLambda$fiUser),
  Y = constants$earthRadius * cos(userFiLambda$lambdaUser) * sin(userFiLambda$fiUser),
  Z = constants$earthRadius * sin(userFiLambda$lambdaUser)
)
print(userCoords)
##          X        Y        Z
## 1 4340.767 4391.267 1597.936

Alpha and beta

# Klobuchar ionospheric model constants
ionisation <- data.frame(
  alphaIon <- c(1.4900e-08, -7.4510e-09, -5.9600e-08, 1.1920e-07),
  betaIon  <- c(1.2900e+05, -1.9660e+05, 6.5540e+04, 3.2770e+05)
)
  
print(ionisation$alphaIon)
## [1]  1.490e-08 -7.451e-09 -5.960e-08  1.192e-07
print(ionisation$betaIon)
## [1]  129000 -196600   65540  327700

sp3 data

# Defining variables used in Klobuchar model and loading data
firstDate <- NA
sp3DataFormatted <- list()
sp3Data <- read.delim('data\\data.SP3', header=FALSE)

# Looping over data (skipping first 23 lines which are header data) and loading it into an array
for (rowIndex in 23:nrow(sp3Data)-1) {
  # Loading data from sp3Data at index 'rowIndex' and splitting the contents (after removing whitespaces) into an array
  # gsub replaces all types of whitespaces with ' '
  # str_trim removes leading and trailing whitespaces
  # str_split splits the content into an array depending on ' '
  row <- str_split(str_trim(gsub('\\s+', ' ', sp3Data[rowIndex, ])), ' ')[[1]]

  # If row length is 7 that line contains timestamp for the next 78 observations
  # Timestamp format: * year month day hour minutes seconds
  if (length(row) == 7) {  
    # Load timestamp into a variable
    dateTime <- strptime(paste(row[2], row[3], row[4], row[5], row[6]), format = "%Y %m %d %H %M")
    
    # Save first observation into a separate variable for later calculations
    if (is.na(firstDate)) firstDate <- dateTime
    
    # gpsTime variable contains difference between first reading and current reading
    gpsTime <- as.numeric(difftime(dateTime, firstDate, units='mins'))
    
    # Skips to the next step of the for loop
    next
  }
  
  # All other rows contain observations for each satellite
  # Observation format: id x y z clock deviations
  data = data.frame(DateTime=dateTime,
                    GPSTime=gpsTime,
                    X=as.numeric(row[2]),
                    Y=as.numeric(row[3]),
                    Z=as.numeric(row[4]))
  
  # Save observation into an array
  sp3DataFormatted[[row[1]]] <- rbind(sp3DataFormatted[[row[1]]], data)
}

# Cleaning unneeded variables from the current workspace
rm(rowIndex, row, data, dateTime, gpsTime, firstDate)

Data calculations and saving

# Calculate Klobuchar ionospheric delay
getIonDist <- function(gpsData){
  # Difference between the user and satellite height
  elevation <- gpsData$Z - userCoords$Z
  
  # Difference between the satellite and user coordinates
  # atan() function is used to compute arctangent angle
  azimuth <- atan((gpsData$X - userCoords$X) - (gpsData$Y - userCoords$Y))
  
  # Angle between the Earth and satellite
  psi <- pi / 2 - elevation - asin(constants$earthRadius / (constants$earthRadius + constants$ippHeight) * cos(elevation))
  # Determine the latitude of the IPP point
  fiIon <- asin(sin(userFiLambda$fiUser) * cos(psi) + cos(userFiLambda$fiUser) * sin(psi) * cos(azimuth))
  # Determine the longitude of the IPP point
  lambdaIon <- userFiLambda$lambdaUser + psi * sin(azimuth) / cos(fiIon)
  # Determine the geomagnetic latitude of the IPP point
  fiMag <- asin(sin(fiIon) * sin(constants$fiPole) + cos(fiIon) * cos(constants$fiPole) * cos(lambdaIon - constants$lambdaPole))
  
  # Determine the local time in the IPP point
  time <- 43200 * lambdaIon / pi + gpsData$GPSTime
  if (time >= 86400) time <- time - 86400
  if (time < 0) time <- time + 86400
  
  # Calculate the amplitude of the ionospheric delay
  azimuthIon <- 0
  for (i in seq(1, 4, 1)) azimuthIon <- azimuthIon + ionisation$alphaIon[i] * `^`(fiMag / pi, i - 1)
  if (azimuthIon < 0) azimuthIon <- 0
  
  # Calculate the period of the ionospheric delay
  psiIon <- 0
  for (i in seq(1, 4, 1)) psiIon <- psiIon + as.double(ionisation$betaIon[i]) * `^`(fiMag / pi, i - 1)
  if (psiIon > 72000) psiIon <- 72000
  
  # Calculate the phase of the ionospheric delay
  xIon <- 2 * pi * (time - 50400) / psiIon
  
  # Calculate ionospheric mapping function (inclination factor)
  Fun <- `^`(1 - `^`(constants$earthRadius / (constants$earthRadius + constants$ippHeight) * cos(elevation), 2), -1 / 2)
  
  # Calculate the value of vertical ionospheric delay
  if (abs(xIon) < pi / 2) dIon <- (5e-9 + azimuthIon * cos(xIon)) * Fun
  else dIon <- 5e-9 * Fun
  
  # Calculate the distance from ionospheric delay by multiplying ionospheric delay with the speed of light
  dIon <- dIon * constants$lightSpeed
  
  # Calculate the Euclidean distance in km from the actual distance between the user and the GPS satellite
  diffXSquared <- (userCoords$X - gpsData$X) ** 2
  diffYSquared <- (userCoords$Y - gpsData$Y) ** 2
  diffZSquared <- (userCoords$Z - gpsData$Z) ** 2
  euclidianDistance <- sqrt(diffXSquared + diffYSquared + diffZSquared)
  
  # Calculate the pseudo distance(m) from the Euclidean distance(km) and ionospheric distance(m)
  distance <- euclidianDistance * 1000 + dIon
  
  return(list(dIon=dIon, distance=distance))
}
# Function for saving and plotting data
saveData <- function(gpsName, lines, ionVector, distVector){
  time <- seq(0, 24, len=length(ionVector))
  
  # Create files for the current GPS
  fileIon  <- file(paste('results/', gpsName, 'ionOutput.txt'))
  fileDist <- file(paste('results/', gpsName, 'distOutput.txt'))
  
  # Save data into the file
  writeLines(as.character(ionVector), fileIon)
  close(fileIon)
  writeLines(lines, fileDist)
  close(fileDist)
  
  # Plotting data 
  plot(time, ionVector, type='l', main=str_interp(paste('Klobucharov model za', gpsName)), xlab='time[h]', ylab='dIon[m]')
  plot(time, distVector, type='l', main=str_interp(paste('Pseudoudaljenost za', gpsName)), xlab='time[h]', ylab='distance[m]')
  
  # Saving plots
  png(paste('plots/', gpsName, 'ionPlot.png'), width=20, height=10, units='cm', res=500)
  plot(time, ionVector, type='l', main=str_interp(paste('Klobucharov model za', gpsName)), xlab='time[h]', ylab='dIon[m]')
  dev.off()
  png(paste('plots/', gpsName, 'distPlot.png'), width=20, height=10, units='cm', res=500)
  plot(time, distVector, type='l', main=str_interp(paste('Pseudoudaljenost za', gpsName)), xlab='time[h]', ylab='distance[m]')
  dev.off()
}
# Preparing data to save
emptySpaces     <- strrep(' ', 25)
timestampHeader <- substr(str_interp('TIME [min]${emptySpaces}'), 0, 15)
distanceHeader  <- substr(str_interp('DISTANCE [m]${emptySpaces}'), 0, 15)
idHeader        <- substr(str_interp('ID ${emptySpaces}'), 0, 10)
header          <- str_interp('${timestampHeader} ${idHeader} ${distanceHeader}')

# Looping over each satellite
for(gpsName in names(sp3DataFormatted)){
  lines <- header
  ionVector  <- c()
  distVector <- c()
  
  # Looping over data for each satellite and formatting the output
  for(gpsIndex in 1:nrow(sp3DataFormatted[[gpsName]])){
    gpsData <- sp3DataFormatted[[gpsName]][gpsIndex, ]
    
    ionDist    <- getIonDist(gpsData)
    ionVector  <- append(ionVector, ionDist$dIon)
    distVector <- append(distVector, ionDist$distance)
    
    gpsTimeFormatted  <- substr(str_interp('${gpsData$GPSTime}${emptySpaces}'), 0, 15)
    distanceFormatted <- substr(str_interp('${ionDist$distance}${emptySpaces}'), 0, 15)
    gpsNameFormatted  <- substr(str_interp('${gpsName}${emptySpaces}'), 0, 10)
    lines <- append(lines, str_interp('${gpsTimeFormatted} ${gpsNameFormatted} ${distanceFormatted}'))
  }
  
  # Saving data using the previously defined function
  saveData(gpsName, lines, ionVector, distVector)
  print(paste('Finished', gpsName))
}

## [1] "Finished PG01"

## [1] "Finished PG02"

## [1] "Finished PG03"

## [1] "Finished PG04"

## [1] "Finished PG05"

## [1] "Finished PG06"

## [1] "Finished PG07"

## [1] "Finished PG08"

## [1] "Finished PG09"

## [1] "Finished PG10"

## [1] "Finished PG11"

## [1] "Finished PG12"

## [1] "Finished PG13"

## [1] "Finished PG14"

## [1] "Finished PG15"

## [1] "Finished PG16"

## [1] "Finished PG17"

## [1] "Finished PG18"

## [1] "Finished PG19"

## [1] "Finished PG20"

## [1] "Finished PG21"

## [1] "Finished PG22"

## [1] "Finished PG23"

## [1] "Finished PG24"

## [1] "Finished PG25"

## [1] "Finished PG26"

## [1] "Finished PG27"

## [1] "Finished PG28"

## [1] "Finished PG29"

## [1] "Finished PG30"

## [1] "Finished PG31"

## [1] "Finished PG32"

## [1] "Finished PR01"

## [1] "Finished PR02"

## [1] "Finished PR03"

## [1] "Finished PR04"

## [1] "Finished PR05"

## [1] "Finished PR07"

## [1] "Finished PR08"

## [1] "Finished PR09"

## [1] "Finished PR11"

## [1] "Finished PR12"

## [1] "Finished PR13"

## [1] "Finished PR14"

## [1] "Finished PR15"

## [1] "Finished PR16"

## [1] "Finished PR17"

## [1] "Finished PR18"

## [1] "Finished PR19"

## [1] "Finished PR20"

## [1] "Finished PR21"

## [1] "Finished PR24"

## [1] "Finished PR25"

## [1] "Finished PE02"

## [1] "Finished PE03"

## [1] "Finished PE04"

## [1] "Finished PE05"

## [1] "Finished PE07"

## [1] "Finished PE08"

## [1] "Finished PE09"

## [1] "Finished PE10"

## [1] "Finished PE11"

## [1] "Finished PE12"

## [1] "Finished PE13"

## [1] "Finished PE14"

## [1] "Finished PE15"

## [1] "Finished PE18"

## [1] "Finished PE19"

## [1] "Finished PE21"

## [1] "Finished PE24"

## [1] "Finished PE25"

## [1] "Finished PE26"

## [1] "Finished PE27"

## [1] "Finished PE30"

## [1] "Finished PE31"

## [1] "Finished PE33"

## [1] "Finished PE34"

## [1] "Finished PE36"
# Cleaning variables from the current workspace
rm(distanceFormatted, distanceHeader, emptySpaces, gpsIndex, gpsName, gpsTimeFormatted, gpsNameFormatted, idHeader, header, ionVector, distVector, lines, timestampHeader)